// LennardJones.cpp: archivo de proyecto principal.

#include <complex>
#include <fstream>
#include <iostream>
#include <string.h>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include<time.h>


using namespace std;

#define npart 20
#define kB 8.314e-03
#define pi 3.1415926535898
#define avogadro 602.205

template<typename T>inline int anint(T x)
	{  return(int(floor(x+0.5)));  }


#include "ewaldcls.h"

int main(void)
{
	void loadconfig(double **pos, double *lbox);
	double **pos;
	double lbox;
	double kappa=0.001;
	double edirect;

	pos = new double*[npart];
	for (int i=0;i<npart;i++)
	{
		pos[i]=new double[3];
	}
	loadconfig(pos, &lbox);

	int puntos=100;
	double desde=0.5;
	double hasta=4.0;
	double paso = (hasta-desde)/double(puntos);
	double ereal;


	cout.precision(4);
	cout<<"kappa"<<setw(20)<<"Directa"<<setw(20)<<"Reciproca"<<setw(20)<<"Total"<<endl;
	int i=0;
	//for (int i=0;i<puntos;i++)
	//{
		kappa=desde+i*paso;
		kappa=8.0;
		ewald *ewd = new ewald(kappa,lbox, pos);
		//ewd->setup();
		ewd->setupdip();
		ewd->r_ewaldDip();
		ewd->k_ewaldDip();
		//ewd->r_ewaldvar();
		//ewd->r_ewald();
		//ewd->k_ewald();

		edirect=ewd->edirectaDipolar();
		ereal=ewd->edirectaDipolarReplicas();
		cout<<"kappa: "<<kappa<<endl;
		cout<<"self: "<<ewd->selfk<<endl;
		cout<<"Vreal: "<<ewd->vr<<"  Vrecip: "<<ewd->vk<<"  Vtotal: "<<ewd->vr+ewd->vk+ewd->selfk<<endl;
		cout<<"Edirecta: "<<edirect<<endl;
		cout<<"Ereal: "<<ereal<<endl;
		//cout<<i<<setw(10)<<kappa<<setw(15)<<ewd->vr<<setw(15)<<ewd->vk<<setw(15)<<ewd->vr+ewd->vk<<setw(15)<<endl;
		delete ewd;
	//}

	//cout<<"Parte directa: "<<ewd->vr<<endl;
	//cout<<"Parte reciproca: "<<ewd->vk<<endl;
	//cout<<"Energia total: "<<ewd->vr+ewd->vk<<endl;




	for (int i=0;i<npart;i++)
	{
		delete pos[i];
	}
	delete []pos;

}

void loadconfig2(double **pos, double *lbox)
{
	int np,rest, ncorr;
	double lb,phi,rcut,Dt,trelaj,time,paso;
	char dummy[150];
	double vel[3];

cout<<"sss"<<endl;
	ifstream MiArchivo ("/home/jonk/projects/ewaldDipolarLast/config.txt",ios::in);
	char line[150];
	//char dummy[20];
	if(MiArchivo.fail())
	{
		cout<<"No se pudo encontrar el archivo de posiciones "<<endl;
		exit(1);
	}
	if (MiArchivo.is_open())
	{
		MiArchivo.getline(dummy,150);
		MiArchivo.getline(dummy,150);
		MiArchivo.getline(dummy,150);
		MiArchivo.getline(dummy,150);
		MiArchivo>>np>>phi>>rcut>>Dt>>rest>>ncorr>>trelaj>>lb>>lb>>lb>>time>>paso;
		MiArchivo.getline(dummy,150);
		MiArchivo.getline(dummy,150);
		MiArchivo.getline(dummy,150);
		MiArchivo.getline(dummy,150);
		MiArchivo.getline(dummy,150);
		MiArchivo.getline(dummy,150);
		for(int i=0;i<np;i++)
		{
			for(int k=0;k<3;k++)
			{
				MiArchivo>>pos[i][k];
				cout<<pos[i][k];
			}
			for(int k=0;k<3;k++)
			{
				MiArchivo>>vel[k];
			}
		}
		MiArchivo.close();
	}
	*lbox=lb;
}

void loadconfig(double **pos, double *lbox)
{
	double lb=1.0;

	for(int i=0;i<npart;i++)
	{
		for(int k=0;k<3;k++)
		{
			pos[i][k]=double(rand())/RAND_MAX-0.5;
		}
	}

	*lbox=lb;
}
